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SUMMARY 

A specialized, microstructural latti ce model , termed MCFE T for combined 
Monte Carlo Finite Element Technique, has been developed which simulates micro- 
structural evolution in material systems where modulated phases occur and the 
directionality of the modulation is influenced by internal and external stresses. 
In this approach the microstructure is discretized onto a fine lattice. Each 
element in the lattice is labelled in accordance with its microstructural iden- 
tity. Diffusion of material at elevated temperatures is simulated by allowing 
exchanges of neighboring elements if the exchange lowers the total energy of 
the system. A Monte Carlo approach is used to select the exchange site while 
the change in energy associated with stress fields is computed using a finite 
element technique. 

The MCFET analysis has been validated by comparing this approach with a 
closed-form, analytical method for stress-assisted, shape changes of a single 
particle in an inifite matrix. Sample MCFET analyses for multiparticle prob- 
lems have also been run and in general the resulting microstructural changes 
associated with the application of an external stress are similar to that 
observed in Ni-Al-Cr alloys at elevated temperatures. 


INTRODUCTION 

Since many of the physical properties of materials are determined by 
microstructure, it is important to be able to predict and control microstruc- 
tural development. Recently a microstructural lattice model (refs. 1 and 2) 
has been developed which can in principle incorporate all relevant driving 
forces and kinetic considerations in a single simulation. Unlike molecular 
dynamics, this approach was developed specifically to predict macroscopic 
behavior, not atomistic behavior. 

The essential elements of this approach consist of mapping the continuum 
microstructure onto a discrete lattice and defining interactions and dynamics 
for the lattice points which are analogous to those in a continuous system. 
Each lattice site is assigned a label which corresponds to the microstructural 
identity of that site. The evolution of the microstructure is governed by the 
appropriate dynamics and the relevant driving forces. Nonconserved dynamics 



(ref. 3) applies to systems in which the internal energy is dependent on param- 
eters which are not conserved during the temporal evolution of the system, 
while conserved dynamics (ref. 4) applies when said parameters are conserved. 

In the former instance, dynamical evolution is controlled by a Monte Carlo pro- 
cedure in which a site is randomly selected and its label changed to one of the 
other allowed labels, if the energy of the system is lowered. In the latter 
instance, dynamical evolution is controlled by a Monte Carlo procedure in which 
a pair of neighboring sites is randomly selected and their labels exchanged if 
the energy of the system is lowered. 

In the present paper a further refinement of the latter approach, termed 
MCFET for combined Monte Carlo Finite Element Technique, is developed specifi- 
cally to model materials that contain modulated phases where the directionality 
of the modulation is strongly influenced by the application of an external 
stress, an example being the Ni-Al-Cr alloy system. In this approach a finite 
element technique is used to calculate the strain energy associated with inter- 
nal and external stress fields, while microstructural evolution is controlled 
using a Monte Carlo procedure which employs conserved dynamics. 

The first half of this paper describes the details of the MCFET approach, 
while the second half of the paper is devoted to sample MCFET simulations of 
increasing complexity. 


MCFET APPROACH 

To model microstructural changes via the MCFET approach an n x n two- 
dimensional grid utilizing square elements was employed. For a realistic 
analysis a typical value for n could be as high as 100. The spatial identity 
of each element (lattice site) in the grid is uniquely defined by a global num- 
bering scheme depicted in figure 1. Each element also has a chemical or micro- 
structural identity which serves to define the shape, distribution, and volume 
fraction of all phases present in the microstructure. For a two phase micro- 
structure each element in the grid would be associated with phase P (precipi- 
tate) or phase M (matrix) as shown in figure 1. 

As previously mentioned, microstructural changes were modelled in the 
present analysis using an approach which is often referred to as conserved 
dynamics. In this approach microstructural evolution occurs by exchanges of 
neighboring sites and therefore the total number of elements identified with 
phase P and M in the aforementioned grid remain constant. The location of the 
exchange pair is chosen at random and the success or failure of the exchange 
is governed by the change in energy. An exchange which lowers the energy of 
the system is accepted, otherwise the exchange is rejected. 

As the size of the Monte Carlo grid is finite, the external surfaces or 
boundaries of the grid pose a dilema if the problem of interest is one in which 
bulk behavior dominates, as in the case of an "infinite" solid. To eliminate 
these surface effects, periodic boundary conditions were used in the present 
analysis. In this approach the Monte Carlo grid appears infinite, although 
periodic, as elements on the upper boundary are considered nearest neighbors 
of elements on the lower boundary while elements on the right boundary are con- 
sidered nearest neighbors of elements on the left boundary. 
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The energy criteria, which determines the success of an exchange, will be 
confined to three components in the present analysis; the interfacial or 
surface energy, the mechanical or strain energy of an elastic solid, and the 
thermal energy. The interfacial energy is directly proportional to the total 
surface area (perimeter in this two-dimensional analysis) associated with the 
interface between phases P and M. The mechanical energy arises from strains 
generated in the solid by externally applied loads and internal misfit strains 
when coherent precipitate(s) are present, as in a multiphase microstructure. 

The calculation of the mechanical energy term for an arbitrary microstructure 
encountered in such an analysis is difficult if not impossible to attain via 
classical elasticity, however, such calculations are conceptually quite easily 
attained using finite element techniques, and as previously stated, will be 
used in the present approach. The thermal energy term, kT, is directly propor- 
tional to the temperature of the solid, T. The significance of the thermal 
energy is dependent on the ratio E/kT, where E = Esurface + Estrain is the 
change in internal energy for a given exchange. When E/kT << 0 the probabil- 
ity of a successful exchange is high, but when E/kT >> 0 the probability of a 
successful exchange is low. If E/kT is near zero the exchange probability 
is near 0.5. Computationally, this is attained by accepting an exchange if 
[1 - TANH (E/kT) 3/2 is greater than a random number between 0 and 1.. 

Before leaving the general discussion concerning the calculation of total 
energy, it should be noted that scaling of these three energy components rela- 
tive to one another is often the most difficult part of the problem. In many 
instances, experimental data for surface tension or elastic properties of indi- 
vidual phases are questionable or nonexistent, so one can only approximate or 
assume values for these physical constants. Changes in these parameters can 
profoundly alter the outcome of the MCFET analysis as the dominant energy term 
changes. This point will be examined in more detail later in this paper. 

The concept of a time scale in most lattice models is related to the 
number of exchanges. In this paper one monte carlo time step, MCS, will be 
defined as n squared attempted exchanges where n is the size of the n x n 
monte carlo grid. 


FINITE ELEMENT STRESS ANALYSIS 

Although finite element techniques can be applied to the present problem 
without great difficulty, the enormous number of exchanges to be considered in 
a realistic analysis would require an impractical amount of computer time to 
calculate the mechanical energy term over the entire monte carlo grid for every 
exchange. Therefore one must devise a scheme which limits the computation time 
per exchange. In the present analysis this is accomplished by confining the 
stress analysis to a small region centered about the exchange site and calcu- 
lating the mechanical energy term of this smaller, local grid before and after 
the exchange. The value of such an approach is that while the absolute ener- 
gies of the local grid may be seriously affected by limiting the area over 
which the stress analysis is done, the change in energy associated with the 
exchange is, however, less likely to be flawed, as similar errors are generated 
and then cancel when the difference in the pre- and post-exchange energies is 
calculated. Further, as the change in energy values calculated are only used 
to decide if a given exchange is accepted, the numerical results, and therefore 
the associated error, is not accumulated. The sensitivity of the error as a 
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function of grid size used In the stress analysis will be quantified at a later 
point In the paper. 

The smaller, local monte carlo grid on which the stress analysis Is run Is 
shown In figure 2. It Is generated from the larger, global monte carlo grid 
and Is seen to contain 36 elements with the exchange pair residing In the 4 
central elements. The actual element type used In the finite element analysis 
Is the two-dimensional, simplex triangle, and as shown In figure 2 each monte 
carlo element contains two such simplex triangles. With this element type and 
grid size one must generate and solve a system of linear equations containing 
98 unknowns for each exchange. This can be done In a fraction of a second on 
a mainframe computer and It Is therefore possible to perform a realistic MCFET 
analysis which may require thousands of exchanges. 

The elastic properties, modulus, Poisson's ratio, etc., of each monte 
carlo element are assigned values corresponding to the phase identity for that 
element. To model the misfit strains between precipitate and matrix, a volu- 
metric strain is simulated by replacing the usual thermal strain term in the 
finite element formulation with a numerically equivalent misfit strain. For a 
two phase analysis in which the misfit strain is -0.1 percent, the thermal 
strain term of the matrix elements are set equal to 0.05 percent while the 
thermal strain term of the precipitate elements are set equal to -0.05 percent. 
Note that the thermal strain term referred to here is not related to the afore- 
mentioned thermal energy term. 

The finite element code used in the present analysis is simple but never- 
theless time efficient, and was adopted from the text of Segerlind (ref. 5). 

In this code the nodal displacements of the finite element grid are computed 
by minimizing the potential energy of the system for a given set of boundary 
conditions. From these nodal displacements one then computes the element 
strains and stresses, and finally the mechanical energy of the local monte 
carlo grid. The choice of boundary conditions are somewhat ambiguous in the 
present problem as one is analyzing a small section of an "infinite" solid. 

As previously noted, extreme accuracies of the absolute energies are not neces- 
sary in the present analysis, therefore the simplest and computational fastest 
boundary conditions were adopted. These set of boundary conditions eliminate 
rigid body motion by fixing the displacements of the central node at zero, 
while rigid body rotation is prevented by setting the lateral displacement of 
the node directly above the central node at zero. Otherwise the body is essen- 
tially unconstrained. As previously stated, the exchange pair is always con- 
tained in the four central elements of the local monte carlo grid and as a 
result there is no directional bias to the exchange resulting from the uncon- 
strained or free boundary conditions used herein, this would not be true if the 
exchange pair were not equidistant from all the outer boundaries of the local 
monte carlo grid. In addition to the displacement boundary conditions, if 
externally applied stresses are to be modelled the appropriate loads are 
applied to the nodes along the outer boundaries. For a uniaxial stress state, 
loads of equal but opposite magnitude are applied to the upper and lower bound- 
ary nodes of the finite element grid. 

The actual codes used to perform the combined monte carlo/finite element 
analysis maybe found in the appendix. In addition to these codes there is a 
complete description of the input and output data used in an analysis as well 
as a simplified flowchart describing the overall approach. 
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INTEGRATION OF THE GRIDS 


Although some description concerning integration of the finite element and 
monte carlo grids has preceded this point, it is felt that a separate, more 
detailed explanation is warranted as this process is of central importance to 
the present analysis. The description of this process is most readily handled 
by example. Given the global monte carlo grid of figure 1, one first selects 
a site at random and then selects one of eight adjoining elements to form the 
exchange pair. If these are of different identity the exchange is attempted. 
Taking element E[7,12] in figure 1 as the site for the exchange the next step 
would involve selection of a neighboring element. In the present analysis this 
is done by choosing a number at random between one and eight, inclusive. As 
seen in figure 3, this number defines the exchange pair according to the con- 
vention shown. The four central elements of the 6x6 local monte carlo grid on 
which the stress analysis will, be done are determined by the logic described 
in figure 3. The remainder of the local monte carlo grid is then filled in 
using the information found in the global monte carlo grid. Continuing with 
the present example, suppose an exchange type 6 was selected, the resulting 
grid would appear as shown in figure 2. Note that since the exchange pair is 
located near the right boundary of the global monte carlo grid, figure 1, the 
local monte carlo grid, figure 2, is composed of elements on both the right 
and left hand side of the global monte carlo grid. The finite element based 
stress analysis is then performed on the local grid to yield the pre-exchange 
strain energy. This is then summed with the interfacial energy. The exchange 
is now performed on the local monte carlo grid and the total energy, both 
interfacial and strain, is recomputed for the new configuration. If the dif- 
ference in the energy before and after exchange is less than the thermal energy 
at this position and time, the exchange is accepted and the global monte carlo 
grid is modified by exchanging the identity of sites E C 7 , 1 2 ] and E C 8 , 1 1 ] in 
figure 1. A new site and exchange type is chosen at random and the process is 
repeated. As previously stated, no numeric information concerning energy, 
especially strain energy, is carried forth, only the conf igurational changes 
of the global monte carlo grid affects subsequent calculations of strain and 
interfacial energy terms. 


VALIDATION OF THE MCFET APPROACH 

Before analyzing more complex problems with multiple particle microstruc- 
tures it is instructive to look at some simple one and two particle problems to 
understand and validate the MCFET approach. In the following analyses, thermal 
and surface tension effects will be ignored and the elastic properties of all 
phases will be taken to be isotropic. 

The first and simplest example to be run tracks the path of a single pre- 
cipitate, comprising one monte carlo element, over an extended period of time 
on a small 6x6 grid. The particle has the same elastic properties as the 
matrix and a misfit of -0.1 percent, further there is no applied stress. The 
results of this analysis are summarized in figure 4 after 100 MCS and shows the 
number of exchanges per site involving the single precipitate particle. Two 
points are to be made. First, of the approximately 3600 attempted exchanges 
only 200 involved the precipitate. This number is consistent with random selec- 
tion of the exchange pair. Second and perhaps more important, there is no 
apparent positional bias introduced by the use of periodic boundaries or the 
finite element stress analysis. 
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In figure 5 a second, single element particle has been added to the 6x6 
grid and after only 16 MCS the particles coalesce. Further, in this time only 
four exchanges involving the precipitate particles were accepted, all of which 
tended to diminish the separation of the two particles. This shows that there 
is an attractive force between these small, single element particles. This 
attractive force is also evident when the elastic properties of matrix and pre- 
cipitate are unequal. Examining the energy change associated with these four 
exchanges, it is apparent that particle coalescence produces the most pro- 
nounced decrease in system energy. By comparison, the energy change of the 
other three exchanges was less than 10 percent of that associated with particle 
coalecsence. The presence of this attractive force was not anticipated, as 
classical elasticity theory predicts misfitting particles in a matrix with 
identical properties are neither attracted nor repelled by one another. The 
origin of this force is apparently related to the finite grid size and boundary 
conditions employed. The discrepancy with classical elasticity does not appear 
to render the technique unusable, as will be shown. Rather, the most signifi- 
cant affect of this attraction results from the large energy decrease accompan- 
ing particle coalescence. This produces an inherent surface tension over and 
above that normally associated with the surface energy term routinely employed 
in these lattice models. 

As previously stated, the size of the finite element grid must in general 
be smaller than the monte carlo grid for computational reasons. In this paper 
the size of the finite element grid will be limited to 6x6. To determine if 
this seriously affects the stress analysis a series of finite element analysis 
were run on the microstructure shown in figure 6. In these analyses, the mis- 
fit is -0.1 percent, the moduli of the particles is 8e 1 0 Pa versus lOelO Pa for 
the matrix, and a stress of 3e8 Pa is applied. The change in strain energy for 
the exchange of sites P and M were then computed for a 6x6, 12x12, and 18x18 
finite element grid. While the absolute energies are much different for the 
three grid sizes, as evident by the pre-exchange strain energy, the change in 
strain energies differ by less than 2 percent. Hence, it appears that the 
choice of a 6x6 finite element grid provides sufficient accuracy while minimiz- 
ing computational requirements. 

To validate the MCFET approach, mi crostructural shape changes of a singu- 
lar particle associated with the application of a uniaxial stress field of 
varying magnitude and direction were examined. This problem has been analyzed 
by others using the approach of Eshelby. In particular, Pineau (ref. 6) has 
presented a thorough analysis of the single-particle problem in which the equi- 
librium particle shape is mapped for a variety of situations. This map has 
been reproduced in figure 7(a). The ratio of precipitate to matrix stiffness 
is plotted on the y-axis, while the ratio of applied stress to misfit, normal- 
ized by the matrix stiffness, is plotted on the x-axis. As seen figure 7(a), 
the application of stress can have three effects; the particle may remain 
spherical, the particle may elongate such that the elongated direction is par- 
allel to the stress axis, or the particle may elongate such that elongated 
direction is perpendicular to the stress axis. Also shown on this map are dis- 
crete symbols representing MCFET runs. In these instances the starting micro- 
structures consisted of a single, square particle, containing four elements, as 
shown in figure 7(b). The finite element grid was 6x6 while the monte carlo 
grid was 12x12, thereby assuring that the particle would be essentially con- 
tained in an infinite matrix. The final particle shape which evolved was a 
singular, rectangular particle, 4x1, with one exception. This particle 
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remained square. These shapes were analogous to that predicted by the Eshelby 
approach. The agreement between the two techniques shows the additional Inter- 
element attraction documented In figure 5 did not seriously influence the out- 
come of the MCFET analysis of the single particle problem. 

At this time it is instructive to extend this analysis to the two multi- 
element, particle problem pictured in figure 8. Here one examines the effect 
of interparticle separation, Rp, on particle shape as a function of applied 
stress. This analysis will be limited to a precipitate moduli and misfit of 
8el0 Pa and -0.1 percent respectively, and a matrix moduli of lOelO Pa. Such a 
choice traverses the Pineau map, figure 7(a), at Ep/Em of 0.8. It is apparent 
that the interparticle separation influences the stress level at which the ini- 
tial configuration becomes unstable. The configuration which evolves is still 
in general agreement with Pineau' s map for the single particle problem. While 
one might suspect adjacent particles to have some effect, it is unclear at this 
point how much of said effect is an artifact of the finite element technique 
similar to that associated with the interelement attraction documented in 
figure 5, and how much is a "real" effect that should be modelled. 


SCALING EFFECTS 

To this point the MCFET analyses have been run without temperature or sur- 
face tension effects. The exchange criteria has been based solely on the 
change in strain energy. As a consequence the absolute dimensions of a parti- 
cle are irrelevent. However, when surface tension effects are included, 
assigning an absolute dimension to a particle is very important as surface 
energy effects scale with the square of the particle size (surface area) while 
strain energy effects scale with the cube of the particle size (volume). In 
the remaining MCFET analyses the effects of surface tension and temperature 
will be included. As a starting point a surface tension of 2e-2 J/m^ will be 
assumed, and each monte carlo element will be assigned a length unit of le-6 m 
on edge. This length scale is appropriate for microstructural analysis. 

To introduce temperature effects in the problem, the thermal energy term, 
kT, is included in the energy balance, as described in the section entitled 
MCFET APPROACH. In essence the magnitude of kT is compared to other energy 
terms. In the present analysis, k is equated to the surface tension. The 
rationale for this approach is, it balances the thermal and surface energy 
terms of an exchange which creates a single interface between precipitate and 
matrix when T = 1 . In other words, the solution temperature of the precipi- 
tate is near unity in the absence of strain effects. To demonstrate this, sim- 
ulations were run with kT equal to le-2, 2e-2, and 3e-2 (T = 0.5, 1.0, and 
1.5). In each case the starting microstructure was a single, square particle, 
comprising four elements on a 12x12 monte carlo grid. The misfit was taken to 
be -0.1 percent, the precipitate moduli was 8el0 Pa versus lOelO Pa for the 
matrix, and there was no externally applied stress. As seen in figure 9, at a 
temperature corresponding to kT=le-2 the particle did not break up nor did it 
change shape. At kT=3e-2 the four element particle broke up into four single, 
element particles after 7 MCS. This temperature is obviously at or above the 
solution temperature of the precipitate for the given parameters. At 
kT = 2e-2 the particle breaks up briefly at 8 MCS but rapidly coalesces. How- 
ever, as shown in figure 9, the aspect ratio is continually changing as there 
is a considerable amount of surface mobility at this temperature. 
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As previously stated the choice of a length scale can affect the surface 
area to volume ratio, and therefore the relative magnitude of the various 
energy terms. To illustrate this effect on mi crostructural evolution the pre- 
ceding problem was rerun using a value of kT = 2e-2 and three length scales; 
le-5, le-6, and le-7 m, all of which are still physically reasonable choices. 

As seen in figure 10, on decreasing the length scale to le-7 m the particle 
becomes more prone to break up due to increased mobility of the constituent 
elements. On increasing the length scale to le-5 m, the particle is seen to be 
stable out to 16 MCS. Further, even surface mobility has been eliminated as 
surface tension effects completely overwhelm thermal and strain energy effects. 
This example points out the importance of choosing an appropriate length scale 
and surface tension number for the problem at hand as such choices can pro- 
foundly alter the relative importance of each energy term. 


STRESS ASSISTED MICROSTRUCTURAL EVOLUTION 

In the preceding sections various aspects of the MCFET approach have been 
examined and validated to some degree. In this section, the effects of an 
applied stress on a single and multiple particle problem will be examined under 
more realistic conditions, surface tension and thermal effects will be included. 
Each of these analyses will be run with a 12x12 monte carlo grid initially con- 
taining 2x2 square particle(s) with a misfit of -0.1 percent. A matrix moduli 
of lOelO Pa and a precipitate moduli of 8el0 Pa will be assumed. Further, a 
surface tension value of 2e-2 J/m 2 and kT of 2e-2 will be used with the 
length scale set at le-6 m per element. 

The effect of a tensile and compressive stress on a single square particle 
will be examined first. These results are presented in figure 11 along with 
that for an unstressed particle. In this figure the overall aspect ratio, Ar, 
being the averaged height to width ratio of all particles is plotted as a func- 
tion of simulated time. This definition was adopted as it is possible for the 
particle to break up. With no applied stress the aspect ratio flip flops about 
unity and it appears that a single particle is the most stable configuration as 
it breaks up only briefly at 8 MCS. Under a tensile stress of 3e8 Pa, the 
aspect ratio climbs above unity near 6 MCS and remains there but for a brief 
time at 8 MCS. At about 12 MCS the particle breaks up but the aspect ratio 
still remains above unity. When a compressive stress of -3e8 Pa is applied, 
the particle first breaks up near 4 MCS and shortly thereafter the aspect ratio 
falls below unity where it remains even when the particles coalesce between 6 
and 7 MCS. Although the kinetics of the shape evolution and particle stability 
are different from that obtained when surface and thermal effects are absent, 
figure 7(a), the aspect ratios obtained here for tension and compression are in 
general agreement with the final equilibrium shapes shown in figure 7(a). 

The effect of applied stress on a mul ti parti cle microstructure is examined 
next. The initial mi crostructure , containing approximately 45 vol 7 « of precip- 
itate arranged as a regular array of square particles, is illustrated in 
figure 12. Also shown are the resulting mi crostuctures after 12 MCS when an 
applied stress of 3e8 or -3e8 Pa is present. As seen here, the compressive 
stress produces a microstructure in which alternate layers of precipitate and 
matrix are aligned perpendicular to the stress axis. When the applied stress 
is tensile, a less regular microstructure evolves in which the alignment of 
alternating layers of precipitate and matrix is parallel to the stress axis. 
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The less regular nature of the microstructure under tension was apparent very 
early In the analysis. That this should occur Is not totally surprising as the 
Plneau map shows this type of structure, N region In figure 7(a), Is stable 
over a rather limited area. The continuous, layered microstructure which 
evolved in the simulations is similar to that observed in real systems, such as 
the Ni-Al-Cr alloys. In these alloys the microstructure, prior to the applica- 
tion of an external stress, consists of a regular array of y' precipitates in a 
continuous y matrix, much like the idealized, starting microstructure shown in 
figure 12. On applying an external stress at elevated temperatures near 1000C, 
these alloys form continuous layers of the y and y' phases. This morphology is 
often referred to as a rafted microstructure. For tensile loads, raft forma- 
tion perpendicular to the stress axis have been observed for some alloy compo- 
sitions (ref. 7), while still other alloy compositions form rafts which are 
parallel to the stress axis (ref. 8). The behavior of a particular alloy is 
believed to be dependent on the misfit and elastic properties of precipitate 
and matrix. On applying a compressive load, the sense of the rafting orienta- 
tion has been observed to reverse, as was the case for the MCFET simulation. 


SUMMARY OF RESULTS 

A specialized, microstructural lattice model, MCFET, has been developed 
which simulates microstructural evolution in material systems where modulated 
phases occur and directionality of the modulation is influenced by internal and 
external stresses. The energy of the stress fields is estimated with a finite 
element technique while microstructural evolution is simulated with a Monte 
Carlo approach in which conserved dynamics is used. 

The MCFET analysis has been validated by comparing this approach with a 
closed-form, analytical method for the single-particle, infinite matrix prob- 
lem. Sample MCFET analyses for a multiparticle problem have also been run and 
in general the resulting microstructures are similar to those observed in the 
Ni-Al-Cr alloy system. 

Although the MCFET approach appears to be a promising analytical tool as 
is, much more development is possible. Some areas for future work include: 

1. Improving the accuracy of the finite element stress analysis by deve- 
loping more accurate boundary conditions, use of complex grid elements, and 
development of a variable grid geometry to more accurately represent the dis- 
tant continuum. 

2. Employing three dimensional analysis as opposed to the current two- 
dimensional analysis. 

3. Develop a dislocation simulator, as misfit dislocations are very impor- 
tant to this class of problem. 
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APPENDIX 

The MCFET analysis comprises four computer programs, MGRID, FGRID, MCFET 
and HPLAY. Printouts of each code may be found at the end of this appendix. 
Each of these codes is written in FORTRAN77 and can be run on an IBM 370 
class, mainframe computer. A flowchart shown below describes how these codes 
work together in a MCFET analysis: 



As indicated above MGRID and FGRID are run first. The inputs to MGRID 

are : 

1) Size of the n x n monte carlo grid [MGR0090] 

2) Location of particle [MGR0120] 

3) Size of particle [MGR0150] 

The information in brackets refers to the line number of the appropriate 
statements in MGRID. MGRID continues to ask for the location and size of 
additional particles until the user responds with 0,0 as input to [MGR0120]. 
The output of MGRID is a file (FILEDEF 10) to be used as input by MCFET. 
FGRID has a singular input: 

l)Size of the n x n finite element grid [FGR0050] 

and also produces a file (FILEDEF 11) to be used as input by MCFET. In this 
paper the value of n in FGRID is generally set at 6. 

MCFET is run next and has the following inputs: 


1) Unit length of monte carlo element 

2 ) Temperature 

3) Surface tension 

4) Random number seed and number of finite element calculations 

5) Applied stress, misfit parameter, and size of monte carlo grid 

6) Elastic parameters of matrix D1,D2,D3 

7) Elastic parameters of precipitate particles D1,D2,D3 


[MCF0270 ] 
[MCF0290 ] 
[MCF0310 ] 
[MCF0330 ] 
[MCF0350 ] 
[MCF0370 ] 
[MCF0390 ] 
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Most of the above have been discussed in the body of this paper, however, 
some explanation is in order. First, the dimensions of all inputs must be 
consistent. Second, the random number seed [MCF0330] is used by the random 
number generator GGUBS , which is a predefined IMSL function; the number of 
finite element calculations [MCF0330] is used to limit the computer time per 
run; and the misfit parameter [MCF0350] is equal to -d/2 where d is the alloy 
misfit. Finally, the elastic parameters are defined as follows: 

Dl=Elastic modulus 

D2=0.3*D1 

D3=0.4*D1 

for the analysis employed here. 

After MCFET completes the specified number of finite element 
calculations the resulting microstructure is displayed and the code then 
requests the following additional inputs: 

1) Change the number of finite element calculations [MCF3500] 

2) Number of finite element calculations [MCF3600] 

I 

If one responds with 0 to [MCF3500] the run is terminated, otherwise MCFET 
requests input to [MCF3600] and the process is repeated until a responce of 0 
is provided as input to [MCF3500]. On termination all information concerning 
microstructural evolution is contained in two files ( FILEDEF 12 and 13). 

As indicated in the flowchart HPLAY is then run to view microstructural 
evolution. The inputs to HPLAY are a file produced by MCFET (FILEDEF 12) and 
the following: 

1) Number of finite element calculations and size of the monte [HPL0050] 

carlo grid 

The first item is equal to the last value of said parameter used in MCFET. 



C MONTE CARLO GRID (MGRID) 

INTEGER MGRID(25,25) 

INTEGER R,C 
DO 10 R=1 , 25 
DO 10 C=1 , 25 
MGRID(R,C)=1 
10 CONTINUE 

PRINT *, 'SIZE OF MONTE CARLO GRID' 

READ *,NGS 

PRINT *, 'ORIGIN TAKEN AS LOWER LEFT CORNER' 
20 PRINT * , ' ROW , COLUMN OF PARTICLE ORIGIN' 
READ *,R,C 

IF (R .EQ. 0) GO TO 35 

PRINT *, 'PARTICLE HEIGHT, WIDTH' 

READ *,H,W 
DO 30 I=R,R+H-1 
DO 30 J=C,C+W-1 
MGRID( I , J)=2 
30 CONTINUE 

GO TO 20 
35 CONTINUE 

DO 40 1=1, NGS 
DO 40 J=1,NGS 

WRITE (10,50) I , J,MGRID( I , J) 

50 FORMAT (314) 

40 CONTINUE 

DO 70 R=NGS ,1,-1 

WRITE (6,100) (MGRID(R.C) ,C=1,NGS) 

70 CONTINUE 

100 FORMAT (2413) 

STOP 

END 


[MGR0090] 

[MGR0120] 


f MGRO 150] 


I 
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C GRID GENERATION CODE (FGRID) 

INTEGER Nl( 1250), N2( 1250), N3( 1250), P( 12), Q( 12), NS(6, 1250), LS 
REAL B(3, 6, 1250) , AR2( 1250) 

PRINT INPUT LATTICE SIZE AS N X N* 

READ *, LS [FGR0050] 

NE=( 2*LS)*LS 

NP=(LS+1)**2 

NC=LS 

NR=2*LS 

N4=INT(NP**0. 5) 

N5=2*N4 

P(l)=l 

P( 2)=2 

P( 3)=N4+1 

P(4)=2 

P(5)=N4+2 

P(6)=N4+1 

P( 7 )=2 

P(8)=N4+3 

P(9 )=N4+2 

P(10)=2 

P( 11)=3 

P( 12)=N4+3 

Q( 1 )=N4+1 

Q(2)=N5+2 

Q(3)=N5+1 

Q(4)=N4+1 

Q(5)=N4+2 

Q(6)=N5+2 

Q ( 7 )=N4+2 

Q(8)=N4+3 

Q(9)=N5+2 

Q(10)=N4+3 

Q( 1 1 )=N5+3 

Q(12)=N5+2 

ITERMl=NC/2 

ITERM2-NR/4 

DO 3 1=1 , ITERM1 

DO 3 J=1 , ITERM2 

DO 3 K=l,4 

I l=K+4*( J- 1 )+2*NR*( I - 1 ) 

I2=3*(K-1) 

I 3=2* (I - 1 )*INT(NP**0 . 5 )+2*( J- 1 ) 

IDUM1=1+I2 
IDUM2=2+I2 
IDUM3=3+I2 
N1(I1)=P(IDUM1)+I3 
N2( I 1 )=P( IDUM2)+I3 
N3( I 1 )=P( IDUM3)+I3 
I1=K+4*(J-1)+(2*(I-1)+1)*NR 
N 1 ( 1 1 ) =Q ( IDUM1)+I3 
N2( I 1 )=Q( IDUM2)+I3 
N3( 1 1 )=Q( IDUM3)+I3 
3 CONTINUE 
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DO 4 1=1, NE 
IY1=(N1(I)-1)/(LS+1) 

X1=(N1(I)-IY1*(LS+1))-1 
IY2=(N2( I ) -1)/(LS+1) 

X2=(N2( I ) - IY2*(LS4-1 ) ) -1 
IY3=(N3(I)-1)/ (LS+1) 

X3=(N3( I ) -IY3*(LS+1) ) -1 

Y1=IY1 

Y2=IY2 

Y3=IY3 

B(l, 1, I)=Y2-Y3 
B( 1 , 3 , I)=Y3-Y1 
B( 1 , 5 , I )=Y1-Y2 
B(2,2, I)=X3-X2 
B(2 ,4, I )=X1 -X3 
B(2 , 6 , I )=X2-X1 
B(3, 1 , I )=B(2,2, I) 

B( 3 , 2 , I )=B( 1,1,1) 

B(3, 3, I)=B(2,4, I) 

B(3 ,4, I )=B( 1,3,1) 

B(3,5,I)=B(2,6,I) 

B(3 , 6 , I )=B( 1,5,1) 

AR2(I)=X2*Y3+X3*Y1+X1*Y2-X2*Y1-X3*Y2-X1*Y3 

NS( 1 , I )=N1 ( I )*2- 1 

NS( 2 , I )=N1 ( I )*2 

NS(3, I)=N2(I)*2-1 

NS(4, I)=N2(I)*2 

NS(5,I)=N3(I)*2-1 

NS(6, I)=N3(I)*2 

WRITE (11,6) I , AR2 ( I ) 

WRITE (11,7) NS (1, I), NS (2, I), NS (3, I), NS (4, I), NS(5, I), NS (6, I) 
DO 9 J=1 , 3 
DO 9 K=1 , 6 

WRITE (11,8) B( J,K, I ) 

6 FORMAT (I4,F10.4) 

7 FORMAT (614) 

8 FORMAT (E12.6) 

9 CONTINUE 

4 CONTINUE 

STOP 

END 
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C COMBINED MONTE CARLO FINITE ELEMENT STRESS ANALYSIS (MCFET) 

C THIS VERSION OF THE CODE INCLUDES RANDOM THERMAL EFFECTS 
C THIS VERSION ALSO INCLUDES SURFACE ENERGY EFFECTS 
C MAXIMUM SIZE OF MONTE CARLO GRID 25X25 (MGRID DATA) 

C STRESS ANALYSIS DONE ON A 6X6 FINITE ELEMENT GRID (FGRID DATA) 

C 98 GLOBAL DEGREES OF FREEDOM (NP) 

C 72 TRIANGULAR SIMPLEX ELEMENTS (NE) 

C BANDWIDTH=18 (NBW) 

C 

C INPUT DATA 

C 

DOUBLE PRECISION DSEED 

INTEGER R,C,RO,CO,Rl,Cl,RN,CN,RS,CS,E 

INTEGER M(6, 6) ,MGRID(25 , 25) ,NPHAS( 72) ,NS(6, 72) 

REAL GF(98),GU(98) , GSM (9 8, 98) , AR2( 72) ,RN3( 3) 

REAL B(3,6,72),D(3,3,2),CM(6,3) ,ESM(6 , 6) ,EF(6) , STRA( 3) , STRE(3) 

TIM=0 . 0 

T=1 . 0 

NP=98 

NE=72 

NBW=18 

NMCS=0 

CALL ERRSET (207,0,-1) 

CALL ERRSET (208,0,-1) 

CALL ERRSET (209,0,-1) 

PRINT *, 'UNIT LENGTH' 

READ *,UL [MCF0270] 

PRINT *, 'TEMPERATURE (0=C0LD TO 1=H0T) ' 

READ *,TEMP [MCF0290] 

PRINT *,' SURFACE ENERGY OF INTERFACE' 

READ * , SURFE [MCF0310] 

PRINT *, 'DSEED, ITERATIONS' 

READ *, DSEED, ITER [MCF03301 

PRINT *,' YLOAD, EPHAS1, SIZE OF MONTE CARLO GRID' 

READ *, YLOAD, EPHASl.NLS 
PRINT *,'D1, D2 , D3 FOR PHASE l’ 

READ *,D1P1,D2P1,D3P1 
PRINT *,'D1, D2, D3 FOR PHASE 2' 

READ *,D1P2,D2P2,D3P2 
EPHAS2=-EPHAS1 
WRITE (12,62) YLOAD, DSEED 
WRITE (12,63) D1P1 ,D2P1 ,D3P1 ,EPHAS1 
WRITE (12,63) D1P2 , D2P2 , D3P2 , EPHAS2 
WRITE (12,64) TEMP, SURFE, UL 

62 FORMAT (2E12.6) 

63 FORMAT (3E12. 6,E12. 6) 

64 FORMAT (3E12.6) 

YLOAD=YLOAD*UL 

C INPUT MGRID DATA 

DO 10 R=1 , NLS 
DO 10 C=1 , NLS 

READ (10,15) IDUM, JDUM,MGRID(R,C) 

WRITE (12,15) IDUM, JDUM,MGRID(R,C) 

15 FORMAT (314) 


[MCF0350] 

[MCF0370] 

[MCF0390] 
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o cj 


10 CONTINUE 

C INPUT FGRID DATA 

DO 20 1=1, NE 
READ (11,25) IDUM, AR2( I ) 

READ (11,35) NS (1,1), NS (2, I), NS (3, I), NS (4, I), NS (5, I), NS (6, I) 
DO 30 J=1 , 3 
DO 30 K=1 , 6 
READ (11,45) B(J,K, I) 

25 FORMAT (I4.F10.4) 

35 FORMAT (614) 

45 FORMAT (E12.6) 

30 CONTINUE 

20 CONTINUE 

C CONSTRUCT MATERIALS PROPERTY MATRIX 

D(1,1,1)=D1P1*UL 
D(2, 2, 1)=D1P1*UL 
D(3, 3, 1)=D3P1*UL 
D( 1,2, 1)=D2P1*UL 
D(2, 1, 1)=D2P1*UL 
D( 1 , 3, 1)=0. 0 
D(3, 1, 1)=0. 0 
D(2, 3, 1 )=0 . 0 
D(3, 2, 1)=0 . 0 
D(l, 1,2)=D1P2*UL 
D(2 , 2, 2)=D1P2*UL 
D(3, 3, 2)=D3P2*UL 
D( 1 , 2, 2)=D2P2*UL 
D(2, 1,2)=D2P2*UL 
D( 1 , 3, 2)=0 . 0 
D(3, 1, 2)=0. 0 
D(2, 3, 2)=0 . 0 
D( 3, 2, 2)=0 . 0 

MAIN CONTROL LOOP 
C 

C SITE AND EXCHANGE SELECTION 

11 IF (NMCS .GE. ITER) GO TO 21 
CALL GGUBS(DSEED, 3,RN3) 

R0=INT(RN3(1)*NLS+1) 

C0=INT(RN3(2)*NLS+1) 

E=INT(RN3(3)*8+1) 

TIM=TIM+1 

R1=R0 

C1=C0 

IF (E .EQ. 1 .OR. E .EQ. 5 .OR. E .EQ. 8) C1=C0+1 

IF (E .EQ. 6 .OR. E .EQ. 3 .OR. E .EQ. 7) C1=C0-1 

IF (E .EQ. 6 .OR. E .EQ. 2 .OR. E .EQ. 5) R1=R0+1 

IF (E .EQ. 7 .OR. E .EQ. 4 .OR. E .EQ. 8) R1=R0-1 

IF (R1 .GT. NLS) R1=R1-NLS 
IF (R1 .LT. 1) R1=R1+NLS 
IF (Cl .GT. NLS) C1=C1-NLS 
IF (Cl .LT. 1) C1=C1+NLS 

IF (MGRID(R0,C0) .EQ. MGRID(R1 ,C1) ) GO TO 11 
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NMCS=NMCS+1 

C CONSTRUCT M ARRAY BEFORE EXCHANGE 
IF (E .EQ. 1 .OR. E .EQ. 5) THEN 
RS=3-R0 
CS=3-C0 
END IF 

IF (E .EQ. 2 .OR. E .EQ. 6) THEN 
RS=3-R0 
CS=4-C0 
END IF 

IF (E .EQ. 3 .OR. E .EQ. 7) THEN 
RS=4-R0 
CS=4-C0 
END IF 

IF (E .EQ. 4 .OR. E .EQ. 8) THEN 
RS=4-R0 
CS=3-C0 
END IF 

DO 40 R=6, 1,-1 
DO 40 C=1 , 6 
RN=R-RS 
CN=C-CS 

IF (RN .GT. NLS) RN=RN-NLS 
IF (RN .LT. 1) RN=RN+NLS 
IF (CN .GT. NLS) CN=CN-NLS 
IF (CN .LT. 1) CN=CN+NLS 
M(R, C)=MGRID(RN,CN) 

40 CONTINUE 

DO 50 R=1 , 6 
DO 50 C=1 , 6 
IDUM=6*R+C-6 
NPHAS(2*IDUM)=M(R,C) 

NPHAS( 2*IDUM-1)=M(R, C) 

50 CONTINUE 

C SURFACE AREA BEFORE EXCHANGE 

SURFO=0 
DO 59 1=3,4 
DO 59 J=2 , 4 

IF (M( I , J) .EQ. M(I,J+1)) SURFO=SURFO+ 1 
IF (M( J, I) .EQ. M(J+1,I)) SURFO=SURFO+ 1 
59 CONTINUE 

SURF0=12-SURF0 
PATH=1 
TMSE=0 . 0 
GO TO 31 

C CONSTRUCT M ARRAY AFTER EXCHANGE 

41 TMSEO=TMSE 
PATH=2 
TMSE=0 . 0 

IF (E .EQ. 1) THEN 
MO=M (3,3) 

M(3,3)=M(3,4) 

M( 3 , 4)=M0 
END IF 
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IF (E .EQ. 2) THEN 
MO=M( 3, 4) 

M(3,4)=M(4,4) 

M(4,4)=M0 
END IF 

IF (E .EQ. 3) THEN 
M0=M(4,4) 

M(4,4)=M(4, 3) 

M(4,3)=M0 
END IF 

IF (E .EQ. 4) THEN 
M0=M(4, 3) 

M(4,3)=M(3,3) 

M(3,3)=M0 
END IF 

IF (E .EQ. 5 .OR. E .EQ. 7) THEN 
MO=M(3, 3) 

M(3,3)=M(4,4) 

M(4,4)=MO 
END IF 

IF (E .EQ. 6 .OR. E .EQ. 8) THEN 
MO=M(3,4) 

M(3,4)=M(4,3) 

M(4,3)=MO 
END IF 
DO 60 R=l,6 
DO 60 C=1 , 6 
IDUM=6*R+C-6 
NPHAS(2*IDUM)=M(R,C) 

NPHAS(2*IDUM-1)=M(R,C) 

60 CONTINUE 

C SURFACE AREA AFTER EXCHANGE 

SURFN=0 
DO 69 1=3,4 
DO 69 J=2,4 

IF (M( I , J) .EQ. M(I,J+1)) SURFN=SURFN+ 1 
IF (M( J, I ) .EQ. M(J+1,I)) SURFN=SURFN+ 1 
69 CONTINUE 

SURFN= 1 2 - SURFN 
GO TO 31 

C COMPARE ENERGY BEFORE AND AFTER EXCHANGE 

51 TMSEN=TMSE 

ECHANG=TMSEN -TMSEO+( SURFN -SURFO ) *SURFE 
NEX=0 

PR0BKT=0 . 5* ( 1 . 0 -TANH ( ECHANG/ ( 2 . 0*TEMP ) ) ) 

CALL GGUBS(DSEED, 1 ,RN3) 

IF (RN3( 1) .LE. PROBKT) THEN 
NEX=1 

MO=MGRID(RO ,C0) 

MGRID(R0,C0)=MGRID(R1,C1) 

MGRID(R1 , C1)=M0 
END IF 

WRITE (12,65) NMCS, RO, CO, Rl, Cl, NEX, ECHANG, TIM 
65 FORMAT (I6.5I4.2E12.6) 
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GO TO 11 


C 

C MECHANICAL ENERGY CALCULATION 

C 

C CONSTRUCT STIFFNESS MATRIX AND FORCE MATRIX 

31 DO 70 1=1, NP 

DO 80 J=1 , NP 
GSM( I , J)=0 . 0 
80 CONTINUE 

GF( I )=0. 0 
70 CONTINUE 

DO 90 KK=1 ,NE 

IF (NPHAS(KK) .EQ. 1) THEN 
ET=EPHAS1 
ELSE 

ET=EPHAS2 
END IF 

C CALCULATE CM=BT X D 

DO 100 1=1,6 
DO 100 J=1 , 3 
CM(I , J)=0. 0 
DO 100 K=1 , 3 

CM( I , J)=CM( I , J)+B(K, I ,KK)*D(K, J,NPHAS(KK) ) 
100 CONTINUE 

C MATRIX MULTIPLICATION TO OBTAIN ESM AND EF 

DO 110 1=1,6 
SUM 1=0.0 
DO 120 K=1 ,2 
SUM1=SUM1+CM( I ,K)*ET 
120 CONTINUE 

DO 110 J=1 , 6 
SUM=0 . 0 
DO 130 K=1 , 3 

SUM=SUM+CM(I,K)*B(K,J,KK) 

130 CONTINUE 

ESM( I , J)=SUM*T/ (2. 0*AR2(KK) ) 

EF( I )=SUMl*T/2 . 0 
110 CONTINUE 

C INSERT ELEMENT MATRIX INTO GLOBAL MATRIX 

DO 140 1=1,6 
II=NS( I ,KK) 

GF( 1 1 )=GF( 1 1 )+EF( I) 

DO 140 J=1 , 6 
JJ=NS( J,KK) 

GSM( I I , JJ)=GSM( I I , JJ)+ESM( I , J) 

140 CONTINUE 

90 CONTINUE 

C ENFORCE BOUNDARY CONDITIONS 

DO 260 1=49,51 
I B= I 

IF (I .EQ. 51) IB=63 
BV=0 . 0 

GSMO=GSM(IB, IB) 

DO 270 J=1 ,NP 



GSM( IB , J)=0. 0 
GSM( J, IB)=0. 0 
270 CONTINUE 

GF( IB)=0. 0 
GSM( IB , IB)=GSMO 
260 CONTINUE 
C APPLY BOUNDARY LOADS 

DO 150 1=1,7 
IYL=2*I 
IYU=2*(42+I ) 

GF( IYL)=GF( IYL) -YLOAD 
GF( I YU)=GF( IYU)+YLOAD 
150 CONTINUE 

C TRIANGULARIZE SYSTEM OF LINEAR EQUATIONS 

DO 160 N=1 ,NP- 1 
NFIN=N+NBW- 1 

IF (NFIN .GT. NP) NFIN=NP 
DO 170 I=N+1 ,NFIN 
RM=GSM( I ,N)/GSM(N,N) 

DO 180 J=N+1 ,NFIN 
GSM( I , J)=GSM( I , J) -GSM(N, J)*RM 
180 CONTINUE 

GF( I )=GF( I ) -GF(N)*RM 
170 CONTINUE 
160 CONTINUE 

C SOLVE LINEAR SYSTEM OF EQUATIONS BY BACKWARD SUBSTITUION 
GU(NP)=GF(NP)/GSM(NP,NP) 

DO 190 I=NP-1, 1,-1 
SUM=0 . 0 
NFIN=I+NBW-1 

IF (NFIN .GT. NP) NFIN=NP 
DO 200 J=I+1,NFIN 
SUM=SUM+GSM ( I , J) *GU ( J ) 

200 CONTINUE 

GU(I)=(GF(I) - SUM) /GSM ( 1,1) 

190 CONTINUE 

C CALCULATE STRESS AND STRAIN 
DO 210 N=1 ,NE 
DO 220 1=1,3 
STRA( I )=0 . 0 
DO 220 K=1 , 6 

STRA( I )=STRA( I )+B(I ,K,N)*GU(NS(K,N) )/AR2(N) 

220 CONTINUE 

DO 230 1=1,3 

IF (NPHAS(N) .EQ. 1) THEN 
ET=EPHAS1 
ELSE 

ET=EPHAS2 
END IF 

IF (I .EQ. 3) ET=0. 0 
STRE(I)=0. 0 
DO 230 K=1 , 3 
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230 

C 


240 

210 

C 


250 

C 


C 

21 


55 

280 


520 


510 

500 


STRE ( 1 )=STRE( I )+D( 1 ,K,NPHAS(N) )*(STRA(K)-ET) 
CONTINUE 

CALCULATE STRAIN ENERGY 
STE=0 . 0 
DO 240 1=1,3 

IF (NPHAS(N) .EQ. 1) THEN 
ET=EPHAS1 
ELSE 

ET=EPHAS2 
END IF 

IF (I .EQ. 3) ET=0 . 0 
STE=STE+STRE ( I )*(STRA( I ) -ET) 

CONTINUE 

TMSE=TMSE+STE*AR2 ( N ) 

CONTINUE 

CALCULATE POTENTIAL ENERGY OF APPLIED LOADS 

SUMPU=0 . 0 

DO 250 1=1,7 

IYL=2*I 

IYU=2*(42+I ) 

SUMPU=SUMPU+YLOAD* ( GU ( IYU) -GU( IYL) ) 

CONTINUE 

CALCULATE TOTAL MECHANICAL ENERGY AND RETURN TO 

TMSE=TMSE/4-SUMPU 

IF (PATH .EQ. 1) GO TO 41 

IF (PATH .EQ. 2) GO TO 51 

OUTPUT FINAL GRID TO SCREEN AND FILE 

PRINT *, 'CHANGE ITER (1=YES 0=N0) ' 

READ *, IDUM1 

PRINT *, ' GRID @ NMCS=',NMCS 
PRINT *, 'TIME=' ,TIM/(NLS*NLS) 

DO 280 R=NLS ,1,-1 

WRITE (6,55) (MGRID(R,C),C=1,NLS) 

FORMAT (2413) 

CONTINUE 

IF ( IDUM1 .EQ. 0) GO TO 520 
PRINT *, ' ITER=' , ITER 
PRINT *, ' NEW VALUE FOR ITER' 

READ *, ITER 
GO TO 11 
DO 500 1=1, NLS 
DO 500 J=1 ,NLS 

WRITE (13,510) I , J,MGRID( I , J) 

FORMAT (314) 

CONTINUE 

STOP 

END 


CONTROL LOOP 


[MCF3500] 


[MCF3600] 
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C PLAYBACK 

INTEGER MGRID(25,25) 

INTEGER R,C,R0,C0,R1 ,C1 

PRINT *, 'ITERATIONS, SIZE OF GRID' 

READ *, ITER,NGS 

READ (12,62) YLOAD.DSEED 

WRITE (6,62) DSEED, YLOAD 

62 FORMAT (2E12.6) 

READ (12,63) D1P1 ,D2P1 ,D3P1 ,EPHAS1 
WRITE (6,63) D1P1 ,D2P1 ,D3P1,EPHAS1 
READ (12,63) D1P2,D2P2,D3P2,EPHAS2 
WRITE (6,63) D1P2,D2P2,D3P2,EPHAS2 

63 FORMAT (3E12. 6,E12. 6) 

READ (12,64) TEMP , SURFE , UL 
WRITE (6,64) TEMP, SURFE, UL 

64 FORMAT (3E12.6) 

DO 10 R=1 ,NGS 
DO 10 C=1 ,NGS 

READ (12,20) IDUM, JDUM,MGRID(R,C) 

20 FORMAT (314) 

10 CONTINUE 

DO 30 R=NGS, 1,-1 

WRITE (6,40) (MGRID(R,C),C=1,NGS) 

40 FORMAT (2413) 

30 CONTINUE 

70 READ (12,60) NMCS ,R0 , CO ,R1 , Cl , NEX,ECHANG,TIM 
60 FORMAT (16,514, 2E12. 6) 

IF (NEX .EQ. 0) GO TO 130 
MO=MGRID(RO,CO) 

MGRID(R0,C0)=MGRID(R1 ,C1) 

MGRID(R1,C1)=M0 

PRINT *, 'NMCS=' ,NMCS, 'E CHANGE=' .ECHANG, 'TIME=' ,TIM 
DO 100 R=NGS ,1,-1 

WRITE (6,110) (MGRID(R,C),C=1,NGS) 

110 FORMAT (2413) 

100 CONTINUE 

130 IF (NMCS .GE. ITER) GO TO 120 

IF (NEX .EQ. 0) GO TO 70 
GO TO 70 
120 STOP 

END 


[IIPL0050] 
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FIGURE 1. - GLOBAL NUMBERING SCHEME USED TO IDENTIFY INDIVIDUAL FIGURE 2. - LOCAL MONTE CARLO GRID GENERATED FROM GLOBAL GRID, FIG- 

MONTE CARLO ELEMENTS. GRID TYPIFIES THAT OF A LOW VOLUME FRAC* URE 1. NOTE CENTRAL LOCATION OF EXCHANGE PAIR P-M, E [7, 12] 

TION, TWO PHASE MICROSTRUCTURE. AND E [8, 11] OF GLOBAL GRID. ALSO SHOWN IS FINITE ELEMENT GRID 

COMPOSED OF 2D SIMPLEX TRIANGLES. 
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FOUR CENTRAL ELEMENTS OF 
LOCAL GRID 


FIGURE 3. - LOGIC USED TO BUILD LOCAL MONTE CARLO GRID FROM GLOBAL 
GRID GIVEN SITE 0 AND AN EXCHANGE TYPE, 1 THRU 8. 
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FIGURE 4. - RANDOM WALK EXPERIMENT FOR A SINGLE ELEMENT PARTICLE. 
THE NUMBER OF EXCHANGES, N, INVOLVING THE PRECIPITATE FOR EACH 
OF THE 36 LATTICE SITES APPEARS TO SHOW THERE IS NO POSITIONAL 


AE/AE 



0 10 20 30 40 50 


TIME (MCS) 

FIGURE 5. - ELEMENT INTERACTION ANALYSIS. NOTE THERE IS AN 
ATTRACTIVE FORCE BETWEEN THE TWO ELEMENTS. FURTHER, ELE- 
MENT COALESCENCE PRODUCES THE GREATEST ENERGY CHANGE, 






FIGURE 6. - VARIATION OF INITIAL STRAIN ENERGIES. E N ', AND CHANGE IN 
STRAIN ENERGIES. AE N ', ASSOCIATED WITH THE INTERCHANGE OF SITES P 
AND II AS A FUNCTION OF THE FINITE ELEMENT GRID SIZE. N. 
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FET ANALYSIS 




FIGURE 7. - COMPARISON OF FINITE ELEMENT ANALYSIS, FET. WITH 
PINEAU'S MAP. BASED ON ESHELBY ANALYSIS, FOR STRESS INDUCED 
SHAPE CHANGES OF A SINGLE PARTICLE, (A). SHAPE EVOLUTION 
OBSERVED IN MCFET SIMULATIONS OF SINGLE PARTICLE ANALYSIS, (B). 
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FIGURE 8. - EFFECT OF INTERPARTICLE SEPARATION, R p , ON PARTICLE SHAPE 
AS A FUNCTION OF APPLIED STRESS, 0. 
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FIGURE S. - SCALING OF THE THERMAL ENERGY TERM. (A) AS kT 
IS INCREASED THERE IS A GREATER TENDENCY FOR THE PART- 
ICLE TO BREAK UP. (B) AT kT = 2x10‘ 2 THE PARTICLE SHAPE 
(H/W) CONTINUALLY CHANGES. 
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FIGURE 10. - CHANGING THE LENGTH UNIT DRAMATICALLY ALTERS 
THE IMPORTANCE OF SURFACE ENERGY BY CHANGING THE SURFACE 
AREA TO VOLUME RATIO. SMALLER PARTICLES (DECREASING UNIT 
LENGTH) HAVE MORE TENDANCE TO BREAK UP, I.E., THEY TEND 
TO DISSOLVE MORE READILY. 





FIGURE 11. - EFFECT OF APPLIED STRESS ON PARTICLE ASPECT 
RATIO, A r . WHEN N p = 1 THE ASPECT RATIO IS THE HEIGHT 
TO WIDTH RATIO OF THE SINGLE PARTICLE. WHEN N p >1, THE 
ASPECT RATIO IS AN AVERAGED HEIGHT TO WIDTH RATIO OF ALL 
PARTICLES. IN THIS ANALYSIS, AN H/W RATIO GREATER THAN 
1 CORRESPONDS TO PARTICLE ELONGATION PARALLEL TO THE 
APPLIED STRESS AXIS. 
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FIGURE 12. - EFFECT OF APPLIED STRESS ON A MULTI-PARTICLE MICROSTRUCTURE. 
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